Now we will use all of the skills we have learned in this workshop in a real world case study. This example walks through using R and publicly available datasets to examine how drought is influencing wind erosion in the Chihuahuan Desert. Our goal is to reproduce Figure 4 from McCord et al. 2023 (https://doi.org/10.1002/ael2.20120). First, we will query the Landscape Data Commons for data that fall within the Chihuahuan Desert ecoregion. Then, we will use these data to set a benchmark of expected wind erosion risk in the Chihuahuan Desert. With this benchmark in place, we will apply it to the monitoring data in the Landscape Data Commons and compare those results at multiple spatial scales. Finally, we will incorporate drought and remote-sensing based fractional cover data, to provide additional context to the field based monitoring data.
Set your working directory make sure it ends in Part 2. If you tab over with a “” Rstudio will help you set it.
getwd() # check your current working drive
## [1] "C:/Users/gharrison/OneDrive - USDA/Documents/R scripts/RinRange_workshop_srm24/RinRangelands-main_1.31.24/RinRangelands-main/Part2"
# setwd("") # change this if needed
Library the packages we will use in this demo
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.3 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.3 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(trex) # for fulling data on the Landscape Data Commons
library(sf) # for processing spatial data
## Linking to GEOS 3.11.2, GDAL 3.7.2, PROJ 9.3.0; sf_use_s2() is TRUE
library(httr2) # for Climate Engine API requests
library(mapview) # for visualizing spatial data
## Warning: package 'mapview' was built under R version 4.3.2
library(ggpubr) # for combining multiple ggplot objects
library(maps)
## Warning: package 'maps' was built under R version 4.3.2
##
## Attaching package: 'maps'
##
## The following object is masked from 'package:purrr':
##
## map
Now we need to pull in some polygons for our region of interest so that we can set a spatial boundaries on our data.
# Bring in the polygons for our region of interest polygons
chihuahuan_desert <- sf::read_sf("MapLayers/ChihuahuanDesert.shp")
sandy_esg <- sf::read_sf("MapLayers/SandyESG.shp")
# extract data from LDC based on the chihuahuan_desert polygon
data <- fetch_ldc_spatial(chihuahuan_desert, data_type = "indicators")
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
head(data)
## rid PrimaryKey
## 1 50790 06080810205047152005-07-19
## 2 50791 06080810212461432005-07-19
## 3 50792 06080810220026862005-07-19
## 4 50793 06080810224129002005-07-20
## 5 50794 06051514053759212006-04-19
## 6 50795 06051514053759212007-03-09
## DBKey ProjectKey
## 1 MURV_2003topresent_DIMA5pt5_11May2021ANON_DatesEdited MURV
## 2 MURV_2003topresent_DIMA5pt5_11May2021ANON_DatesEdited MURV
## 3 MURV_2003topresent_DIMA5pt5_11May2021ANON_DatesEdited MURV
## 4 MURV_2003topresent_DIMA5pt5_11May2021ANON_DatesEdited MURV
## 5 MURV_2003topresent_DIMA5pt5_11May2021ANON_DatesEdited MURV
## 6 MURV_2003topresent_DIMA5pt5_11May2021ANON_DatesEdited MURV
## DateVisited EcologicalSiteId PercentCoveredByEcoSite
## 1 2005-07-19T00:00:00.000Z R042XB010NM NA
## 2 2005-07-19T00:00:00.000Z R042XB035NM NA
## 3 2005-07-19T00:00:00.000Z R042XB010NM NA
## 4 2005-07-20T00:00:00.000Z R042XB027NM NA
## 5 2006-04-19T00:00:00.000Z R042XB021NM NA
## 6 2007-03-09T00:00:00.000Z R042XB021NM NA
## Latitude_NAD83 Longitude_NAD83 LocationStatus LocationType
## 1 33.0389 -107.336 NA NA
## 2 33.0385 -107.337 NA NA
## 3 33.0386 -107.336 NA NA
## 4 32.9359 -107.26 NA NA
## 5 32.563 -106.521 NA NA
## 6 32.563 -106.521 NA NA
## Latitude_NAD83_Actual Longitude_NAD83_Actual BareSoilCover TotalFoliarCover
## 1 NA NA 8.50 33.00
## 2 NA NA 3.00 32.00
## 3 NA NA 7.50 17.00
## 4 NA NA 7.00 57.00
## 5 NA NA 3.75 48.75
## 6 NA NA NA NA
## AH_AnnGrassCover AH_ForbCover AH_GrassCover AH_PerenForbCover
## 1 0.0 23.0 7.5 5.0
## 2 0.0 15.5 12.5 4.5
## 3 0.0 12.0 2.0 7.0
## 4 14.5 8.0 29.5 2.5
## 5 0.0 7.5 2.5 0.0
## 6 NA NA NA NA
## AH_PerenForbGrassCover AH_PerenGrassCover AH_ShrubCover FH_CyanobacteriaCover
## 1 12.5 7.5 5 0
## 2 16.0 12.5 8 0
## 3 8.0 2.0 5 0
## 4 19.0 16.5 17 0
## 5 2.5 2.5 0 0
## 6 NA NA NA 0
## FH_DepSoilCover FH_DuffCover FH_EmbLitterCover FH_HerbLitterCover
## 1 0 0 0 10.5
## 2 0 0 0 8.5
## 3 0 0 0 8.0
## 4 0 0 0 7.5
## 5 0 0 0 17.5
## 6 0 0 0 NA
## FH_LichenCover FH_MossCover FH_RockCover FH_TotalLitterCover
## 1 0 0 48.00 10.50
## 2 0 0 56.50 8.50
## 3 0 0 67.50 8.00
## 4 0 0 28.50 7.50
## 5 0 0 28.75 18.75
## 6 NA 0 NA NA
## FH_VagrLichenCover FH_WaterCover FH_WoodyLitterCover GapCover_101_200
## 1 0 0 0.00 NA
## 2 0 0 0.00 NA
## 3 0 0 0.00 NA
## 4 0 0 0.00 NA
## 5 0 0 1.25 NA
## 6 0 0 NA NA
## GapCover_200_plus GapCover_25_50 GapCover_25_plus GapCover_51_100
## 1 NA NA NA NA
## 2 NA NA NA NA
## 3 NA NA NA NA
## 4 NA NA NA NA
## 5 NA NA NA NA
## 6 NA NA NA NA
## Hgt_Forb_Avg Hgt_Grass_Avg Hgt_Herbaceous_Avg Hgt_PerenForb_Avg
## 1 NA NA NA NA
## 2 NA NA NA NA
## 3 NA NA NA NA
## 4 NA NA NA NA
## 5 NA NA NA NA
## 6 NA NA NA NA
## Hgt_PerenForbGrass_Avg Hgt_PerenGrass_Avg Hgt_Woody_Avg RH_AnnualProd
## 1 NA NA NA <NA>
## 2 NA NA NA <NA>
## 3 NA NA NA <NA>
## 4 NA NA NA <NA>
## 5 NA NA NA <NA>
## 6 NA NA NA <NA>
## RH_BareGround RH_BioticIntegrity RH_CommentsBI RH_CommentsHF RH_CommentsSS
## 1 <NA> <NA> NA NA NA
## 2 <NA> <NA> NA NA NA
## 3 <NA> <NA> NA NA NA
## 4 <NA> <NA> NA NA NA
## 5 <NA> <NA> NA NA NA
## 6 <NA> <NA> NA NA NA
## RH_Compaction RH_DeadDyingPlantParts RH_FuncSructGroup RH_Gullies
## 1 <NA> <NA> <NA> <NA>
## 2 <NA> <NA> <NA> <NA>
## 3 <NA> <NA> <NA> <NA>
## 4 <NA> <NA> <NA> <NA>
## 5 <NA> <NA> <NA> <NA>
## 6 <NA> <NA> <NA> <NA>
## RH_HydrologicFunction RH_InvasivePlants RH_LitterAmount RH_LitterMovement
## 1 <NA> <NA> <NA> <NA>
## 2 <NA> <NA> <NA> <NA>
## 3 <NA> <NA> <NA> <NA>
## 4 <NA> <NA> <NA> <NA>
## 5 <NA> <NA> <NA> <NA>
## 6 <NA> <NA> <NA> <NA>
## RH_PedestalsTerracettes RH_PlantCommunityComp RH_ReprodCapabilityPeren
## 1 <NA> <NA> <NA>
## 2 <NA> <NA> <NA>
## 3 <NA> <NA> <NA>
## 4 <NA> <NA> <NA>
## 5 <NA> <NA> <NA>
## 6 <NA> <NA> <NA>
## RH_Rills RH_SoilSiteStability RH_SoilSurfLossDeg RH_SoilSurfResisErosion
## 1 <NA> <NA> <NA> <NA>
## 2 <NA> <NA> <NA> <NA>
## 3 <NA> <NA> <NA> <NA>
## 4 <NA> <NA> <NA> <NA>
## 5 <NA> <NA> <NA> <NA>
## 6 <NA> <NA> <NA> <NA>
## RH_WaterFlowPatterns RH_WindScouredAreas SoilStability_All
## 1 <NA> <NA> NA
## 2 <NA> <NA> NA
## 3 <NA> <NA> NA
## 4 <NA> <NA> NA
## 5 <NA> <NA> NA
## 6 <NA> <NA> NA
## SoilStability_Protected SoilStability_Unprotected DateLoadedInDb
## 1 NA NA 2023-04-10T00:00:00.000Z
## 2 NA NA 2023-04-10T00:00:00.000Z
## 3 NA NA 2023-04-10T00:00:00.000Z
## 4 NA NA 2023-04-10T00:00:00.000Z
## 5 NA NA 2023-04-10T00:00:00.000Z
## 6 NA NA 2023-04-10T00:00:00.000Z
## mlra_name mlrarsym
## 1 Southern Desertic Basins, Plains, and Mountains 42
## 2 Southern Desertic Basins, Plains, and Mountains 42
## 3 Southern Desertic Basins, Plains, and Mountains 42
## 4 Southern Desertic Basins, Plains, and Mountains 42
## 5 Southern Desertic Basins, Plains, and Mountains 42
## 6 Southern Desertic Basins, Plains, and Mountains 42
## na_l1name na_l2name us_l3name
## 1 NORTH AMERICAN DESERTS WARM DESERTS Chihuahuan Deserts
## 2 NORTH AMERICAN DESERTS WARM DESERTS Chihuahuan Deserts
## 3 NORTH AMERICAN DESERTS WARM DESERTS Chihuahuan Deserts
## 4 NORTH AMERICAN DESERTS WARM DESERTS Chihuahuan Deserts
## 5 NORTH AMERICAN DESERTS WARM DESERTS Chihuahuan Deserts
## 6 NORTH AMERICAN DESERTS WARM DESERTS Chihuahuan Deserts
## us_l4name State modis_landcover horizontal_flux_total_MD
## 1 Chihuahuan Basins and Playas NM Closed Shrublands NA
## 2 Chihuahuan Basins and Playas NM Closed Shrublands NA
## 3 Chihuahuan Basins and Playas NM Closed Shrublands NA
## 4 Low Mountains and Bajadas NM Open Shrublands NA
## 5 Chihuahuan Montane Woodlands NM Grasslands NA
## 6 Chihuahuan Montane Woodlands NM Grasslands NA
## vertical_flux_MD PM2_5_MD PM10_MD
## 1 NA NA NA
## 2 NA NA NA
## 3 NA NA NA
## 4 NA NA NA
## 5 NA NA NA
## 6 NA NA NA
First lets make a map of the region of the boundaries of the chihuahuan desert within the US. To do this, we need to pull in a US state boundary.
# pull in US state boundaries
usa <- st_as_sf(maps::map("state", fill=TRUE, plot =FALSE))
# Generate a map from just the state outlines
state_base_map <- ggplot(usa) +
geom_sf(color = "#2b2b2b", fill = "white", size=0.125)
state_base_map
Good start, lets add the Chihuahuan desert study area on this map. Remember ggplot figures work like layers, so we can add the desert on top of our base map
# lets add in the boundary for the Chihuahan desert
# add this on to our prior map
state_map <- state_base_map +
geom_sf(data = chihuahuan_desert, aes(fill = "Chihuahuan Desert"))+
theme_minimal()
state_map
Nice. Now let’s use some customization within ggplot to adjust legend
position and change colors
us_state_map <- state_base_map +
geom_sf(data = chihuahuan_desert, aes(fill = "Chihuahuan Desert")) + # alpha makes the shape transparent
scale_fill_manual(values = "goldenrod1", guide = guide_legend(title = "Legend")) +
theme(legend.position = "top")
us_state_map
Awesome, now we have a starting map. Lets move on to working with the
data.
Take a look at the LDC data and clean it up
# Make the DateVisited formatted as date
data <- data |> dplyr::mutate(DateVisited = lubridate::as_date(DateVisited),
Year = lubridate::year(DateVisited) |> as.factor())
# Make sure the Lat/Long coordinates are numeric
data <- data |>
dplyr::mutate(Longitude_NAD83 = as.numeric(Longitude_NAD83),
Latitude_NAD83 = as.numeric(Latitude_NAD83))
We want to understand which sites might be susceptible to wind erosion. To do this, we took a subset (~600) of the points in the ecoregion and compared them with wind erosion models and determined that decreases in total foliar cover below 30% are associated with an increase in wind erosion risk. However, we do not want to apply the benchmark to the points used to develop the benchmark. So we will first remove the benchmark generation points. Then we will apply the 30% benchmark to these points.
pull in the benchmark
development data
benchmark <- read.csv("benchmark_points.csv")
head(benchmark)
## PrimaryKey SpeciesState PlotID PlotKey
## 1 1611141556434902016-09-01 NM Bottom-002 1.61114E+14
## 2 16111616394616792016-09-01 NM NR-FIRE-TRT 1611161639461679
## 3 16111616402732582016-09-01 NM NR-FIRE-CTL 1611161640273258
## 4 16102714453525922016-09-01 NM Sandy-463 1610271445352592
## 5 16081110402163562016-09-01 NM Gsand-144 1608111040216356
## 6 16081114554791692016-09-01 NM Hills-188 1608111455479169
## DBKey EcologicalSiteId State County
## 1 NM_LasCrucesDO_LUP_2016_5-3b_01.mdb R042XB012NM NM Dona Ana
## 2 NM_LasCrucesDO_LUP_2016_5-3b_01.mdb R042XB035NM NM Dona Ana
## 3 NM_LasCrucesDO_LUP_2016_5-3b_01.mdb R042XB035NM NM Dona Ana
## 4 NM_LasCrucesDO_LUP_2016_5-3b_01.mdb R042XB011NM NM Luna
## 5 NM_LasCrucesDO_LUP_2016_5-3b_01.mdb R042XB024NM NM Dona Ana
## 6 NM_LasCrucesDO_LUP_2016_5-3b_01.mdb R042XB010NM NM Dona Ana
## DateEstablished DateLoadedInDb ProjectName DateVisited
## 1 11/6/2016 17:00 9/1/2016 New Mexico Las Cruces FO 2016 NA
## 2 11/8/2016 17:00 9/1/2016 New Mexico Las Cruces FO 2016 NA
## 3 11/8/2016 17:00 9/1/2016 New Mexico Las Cruces FO 2016 NA
## 4 10/24/2016 18:00 9/1/2016 New Mexico Las Cruces FO 2016 NA
## 5 8/10/2016 18:00 9/1/2016 New Mexico Las Cruces FO 2016 NA
## 6 8/10/2016 18:00 9/1/2016 New Mexico Las Cruces FO 2016 NA
## source LocationType ELEVATION PercentCoveredByEcoSite MLRA_ID MLRARSYM
## 1 TerrADat <NA> NA NA 43 42B
## 2 TerrADat <NA> NA NA 43 42B
## 3 TerrADat <NA> NA NA 43 42B
## 4 TerrADat <NA> NA NA 43 42B
## 5 TerrADat <NA> NA NA 43 42B
## 6 TerrADat <NA> NA NA 43 42B
## MLRA_NAME LRRSYM LRR_NAME
## 1 Southern Rio Grande Rift D Western Range and Irrigated Region
## 2 Southern Rio Grande Rift D Western Range and Irrigated Region
## 3 Southern Rio Grande Rift D Western Range and Irrigated Region
## 4 Southern Rio Grande Rift D Western Range and Irrigated Region
## 5 Southern Rio Grande Rift D Western Range and Irrigated Region
## 6 Southern Rio Grande Rift D Western Range and Irrigated Region
Start by removing the training benchmark data
# remove data that have been used in benchmarking
nrow(data)
## [1] 2294
data <- data |> subset(!PrimaryKey %in% benchmark$PrimaryKey)
# check total number of plots
nrow(data)
## [1] 1709
Score score the plots based on if they are meeting or not meeting benchmarks
# Apply the benchmark of 30% to the Total Foliar Cover values
data <- data |>
dplyr::mutate(evaluated = dplyr::case_when(
TotalFoliarCover >30 ~ "Meeting",
TotalFoliarCover <=30 ~ "Not meeting")) |>
dplyr::filter(!is.na(TotalFoliarCover))
# this will create a new column called evaluated with values of meeting or not meeting based on the values in the Total Foliar Cover col
# there were some plots with NA in the TotalFoliarCover - we want to remove those since we cannot evaluate them for this benchmark
Convert these data to spatial objects so we can make maps and conduct spatial analyses To accomplish this, we will need specify a CRS (Coordinate reference system) [https://rspatial.org/raster/spatial/6-crs.html]
# move from a st to sf object
# users must specify the columns which are coordinates, and the Coordinate reference system
data_sf <- sf::st_as_sf(data %>% subset(!is.na(Longitude_NAD83)),
coords = c("Longitude_NAD83", "Latitude_NAD83"),
crs = st_crs(4269)) # 4268 is the CRS code for NAD 83
# save these data as a shapefile
sf::st_write(data_sf, "ldc_data.shp", append=F)
## Warning in abbreviate_shapefile_names(obj): Field names abbreviated for ESRI
## Shapefile driver
## Deleting layer `ldc_data' using driver `ESRI Shapefile'
## Writing layer `ldc_data' to data source `ldc_data.shp' using driver `ESRI Shapefile'
## Writing 1502 features with 85 fields and geometry type Point.
We can add points to our base map with states
# recall study area map
us_state_map
This map is great for providing context of our study are in the US…. but
to see the plot points on a map we are going to need to zoom into New
Mexico and AZ. To do that, we can filter the ‘usa’ data frame
head(usa) # seems like the state name is in the ID column
## Simple feature collection with 6 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -124.3834 ymin: 30.24071 xmax: -71.78015 ymax: 42.04937
## Geodetic CRS: +proj=longlat +ellps=clrk66 +no_defs +type=crs
## ID geom
## alabama alabama MULTIPOLYGON (((-87.46201 3...
## arizona arizona MULTIPOLYGON (((-114.6374 3...
## arkansas arkansas MULTIPOLYGON (((-94.05103 3...
## california california MULTIPOLYGON (((-120.006 42...
## colorado colorado MULTIPOLYGON (((-102.0552 4...
## connecticut connecticut MULTIPOLYGON (((-73.49902 4...
zoomed_state_map <- ggplot()+
geom_sf(data=usa %>% dplyr::filter(ID %in% c("new mexico", "arizona")),
fill = "white", color = "black")+
theme_minimal()
zoomed_state_map
Add in the Chihuahan Desert and Sandy ESG map layers to this zoomed in
map
zoomed_state_map <- zoomed_state_map+
geom_sf(data = chihuahuan_desert, aes(fill = "Chihuahuan Desert"), alpha = 0.75) + # make semi transparent with alpha
geom_sf(data = sandy_esg, aes(fill = "Sandy ESG"), color = "gold4" )+
scale_fill_manual(values = c("Chihuahuan Desert" = "goldenrod1",
"Sandy ESG" = "gold4"),
name = "",
labels = c("Chihuahuan Desert", "Sandy ESG")) +
theme(legend.position = "top")
zoomed_state_map
##### challenge!
# add the plot locations to the study area map consider changing dot
size using “size =”
#zoomed_state_map_withPts <- zoomed_state_map +
Now we can see the plot locations, but I wonder where plots are that are meeting or not meeting the benchmark? We can change the color of each plot location based on the ‘evaluated’ column
# get the values in this column
summary(as.factor(data_sf$evaluated))
## Meeting Not meeting
## 948 554
# Specify colors for meeting and not meeting
meeting_colors <- c("Meeting" = "royalblue1", "Not meeting" = "firebrick4")
# add these points to the zoomed state map, but have the color of each point indicate if the location is meeting or not meeting the benchmark
study_with_bench_eval <- zoomed_state_map +
geom_sf(data = data_sf,
size = 0.75,
aes(color = evaluated)) +
scale_color_manual(values = meeting_colors,
name = "Wind Erosion Benchmark")
# Display the map with adjusted point size and color
study_with_bench_eval
Now that we know where these points are, we can make some summary graphs to show the portion of the plots which are meeting this benchmark.
For this, we are interested in understanding how many plots within the Chihuahuan Desert, the Sandy ESG, MLRA 42, and a long-term research site at the Jornada Experimental Range called “NWERN_JER”).
# Identify the membership of each plot within different polygons (Chihuahuan Desert, MLRA 42, SandyESG and a long-term research site at the Jornada Experimental Range called "NWERN_JER")
# We know all the points are within the Chihuahaun desert, do create a col and fill all values with yes
# sort plots by if they occur in MLRA42. If the value in mlrasym col is 42, put 'yes' in the new column called MLRA42
data_sf <- data_sf |>
dplyr::mutate(Chihuahuan_desert = "yes",
MLRA42 = dplyr::case_when(mlrarsym %in% "42"~ "yes"),
NWERN_JER = dplyr::case_when(ProjectKey %in% "NWERN_JER"~ "yes"))
we also want to know if the points falls within the Sandy ESG. We will use the st_within function
# create a col called SandyESG, where values correspond to if each row is within the SandyESG
# Assuming data_sf and sandy_sf are sf objects
data_sf$SandyESG <- sf::st_within(data_sf, sandy_esg)
hmm seems like there is an issue about mismatching CRS
# figure out the current crs for each data
print(st_crs(data_sf))
## Coordinate Reference System:
## User input: EPSG:4269
## wkt:
## GEOGCRS["NAD83",
## DATUM["North American Datum 1983",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Geodesy."],
## AREA["North America - onshore and offshore: Canada - Alberta; British Columbia; Manitoba; New Brunswick; Newfoundland and Labrador; Northwest Territories; Nova Scotia; Nunavut; Ontario; Prince Edward Island; Quebec; Saskatchewan; Yukon. Puerto Rico. United States (USA) - Alabama; Alaska; Arizona; Arkansas; California; Colorado; Connecticut; Delaware; Florida; Georgia; Hawaii; Idaho; Illinois; Indiana; Iowa; Kansas; Kentucky; Louisiana; Maine; Maryland; Massachusetts; Michigan; Minnesota; Mississippi; Missouri; Montana; Nebraska; Nevada; New Hampshire; New Jersey; New Mexico; New York; North Carolina; North Dakota; Ohio; Oklahoma; Oregon; Pennsylvania; Rhode Island; South Carolina; South Dakota; Tennessee; Texas; Utah; Vermont; Virginia; Washington; West Virginia; Wisconsin; Wyoming. US Virgin Islands. British Virgin Islands."],
## BBOX[14.92,167.65,86.45,-40.73]],
## ID["EPSG",4269]]
print(st_crs(sandy_esg))
## Coordinate Reference System:
## User input: WGS 84 / Pseudo-Mercator
## wkt:
## PROJCRS["WGS 84 / Pseudo-Mercator",
## BASEGEOGCRS["WGS 84",
## ENSEMBLE["World Geodetic System 1984 ensemble",
## MEMBER["World Geodetic System 1984 (Transit)"],
## MEMBER["World Geodetic System 1984 (G730)"],
## MEMBER["World Geodetic System 1984 (G873)"],
## MEMBER["World Geodetic System 1984 (G1150)"],
## MEMBER["World Geodetic System 1984 (G1674)"],
## MEMBER["World Geodetic System 1984 (G1762)"],
## MEMBER["World Geodetic System 1984 (G2139)"],
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ENSEMBLEACCURACY[2.0]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]],
## CONVERSION["Popular Visualisation Pseudo-Mercator",
## METHOD["Popular Visualisation Pseudo Mercator",
## ID["EPSG",1024]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["False easting",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["easting (X)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["northing (Y)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["Web mapping and visualisation."],
## AREA["World between 85.06°S and 85.06°N."],
## BBOX[-85.06,-180,85.06,180]],
## ID["EPSG",3857]]
Ah yes the crs for sandy_esg is WGS 84. Let’s change that.
# Transform sandy_esg_transformed to have the same CRS as nm_data
sandy_esg_transformed <- sf::st_make_valid(sf::st_transform(sandy_esg, st_crs(data_sf)))
# try again to assign values based on points, but use the transformed data
data_sf$SandyESG <- sf::st_within(data_sf, sandy_esg_transformed)
head(summary(data_sf$SandyESG), 20) # the values in this col are either 0 or 1
## Length Class Mode
## [1,] 0 -none- numeric
## [2,] 0 -none- numeric
## [3,] 0 -none- numeric
## [4,] 0 -none- numeric
## [5,] 0 -none- numeric
## [6,] 0 -none- numeric
## [7,] 0 -none- numeric
## [8,] 0 -none- numeric
## [9,] 0 -none- numeric
## [10,] 0 -none- numeric
## [11,] 0 -none- numeric
## [12,] 0 -none- numeric
## [13,] 0 -none- numeric
## [14,] 0 -none- numeric
## [15,] 0 -none- numeric
## [16,] 0 -none- numeric
## [17,] 0 -none- numeric
## [18,] 1 -none- numeric
## [19,] 0 -none- numeric
## [20,] 0 -none- numeric
# Convert the logical values to "Yes" or "No" if needed
data_sf$SandyESG <- ifelse(data_sf$SandyESG, "yes", "No")
head(summary(data_sf$SandyESG), 20)
## Length Class Mode
## 1502 character character
Great. Now for each plot (row), we know if the plot occurs within the Chihuahan desert, Sandy ESG, MLRA 42, and in NWERN
To prepare for plotting, we need to transform data from wide to long format using a pivot.
# To faciliate plotting, it would be helpful if the data were long rather than wide.
data_long <- as.data.frame(data_sf) |>
tidyr::pivot_longer(cols = c("MLRA42",
"SandyESG",
"NWERN_JER",
"Chihuahuan_desert"),
names_to = "region",
values_to = "region_membership") |>
subset(region_membership =="yes")|>
dplyr::mutate(region = factor(region)|>
dplyr::recode( "NWERN_JER" = "Site",
"SandyESG" = "ESG",
"MLRA42" = "MLRA",
"Chihuahuan_desert" = "Ecoregion"))
# We also need region to be a factor for easier plot formatting
data_long$region <- factor(data_long$region,
levels = c( "Site",
"ESG",
"MLRA",
"Ecoregion")
)
# box plot
benchmark_boxplot <-
ggplot(data_long |>
subset(Year %in% c("2015", "2016", "2017", "2018", "2019", "2020", "2021", "2022") &
TotalFoliarCover),
aes(x = Year,
y = TotalFoliarCover,
group = Year)) +
facet_grid(rows = "region")+
geom_boxplot()+
geom_point(aes(color = evaluated))+
scale_color_manual(values = c("darkblue", "orange"))+
theme_bw()+
ylab("Total foliar cover (%)")+
guides(color=guide_legend(title="Benchmark status"))
benchmark_boxplot
#benchmark_boxplot_jittered <-
benchmark_barchart <-
ggplot(data_long |>
subset(Year %in% c("2015", "2016", "2017", "2018", "2019", "2020", "2021", "2022") &
TotalFoliarCover),
aes(x = Year,
y = TotalFoliarCover,
fill = evaluated)) +
facet_grid(rows = "region")+
geom_bar(position = "fill", stat = "identity")+
scale_fill_manual(values = c("darkblue", "orange"))+
theme_bw()+
ylab("Proportion of plots")+
guides(fill=guide_legend(title="Benchmark status"))
benchmark_barchart
data_long_summary <-
data_long |>
subset(Year %in% c("2015", "2016", "2017", "2018", "2019", "2020", "2021", "2022") &
TotalFoliarCover) |>
dplyr::select(PrimaryKey, region, Year) |>
dplyr::mutate(Year = as.character(Year)) |>
dplyr::distinct()|>
dplyr::group_by(Year, region) |>
dplyr::tally()
Define Climate Engine API key and run a simple test endpoint to make sure that the key authenticated correctly. The Climate Engine API documentation is published at https://docs.climateengine.org
There are also Climate Engine API tutorials in R and Python here: https://support.climateengine.org/article/42-api-tutorials
The Climate Engine API has a ‘Swagger’ app for testing and making requests at: https://api.climateengine.org/docs
NOTE: The key below will remain active through February 4, please request a free key to run this in the future: https://docs.climateengine.org/docs/build/html/registration.html
# Define root url for Climate Engine API
root_url <- 'https://api.climateengine.org/'
# Define key
key <- 'eyJhbGciOiJIUzI1NiIsInR5cCI6IkpXVCJ9.eyJmcmVzaCI6ZmFsc2UsImlhdCI6MTcwNTg3ODY2NiwianRpIjoiYmEyYTZmMzMtNjdkOS00ODAzLThiOTQtYWQ0Yjg5NmZjZDJkIiwibmJmIjoxNzA1ODc4NjY2LCJ0eXBlIjoiYWNjZXNzIiwic3ViIjoidUZXM0VmMUxhWmI5U0VhUnRPdzVNUzI2UjBTMiIsImV4cCI6MTcwNzE3NDY2Niwicm9sZXMiOiJ1c2VyIiwidXNlcl9pZCI6InVGVzNFZjFMYVpiOVNFYVJ0T3c1TVMyNlIwUzIifQ.WIQZIn97BpFZorL_n6NGvEi_MHg2qMnIUh9Rm6aOcl0'
# Define endpoint
endpoint <- '/home/validate_key'
# Run simple endpoint to get key expiration
test <- request(base_url = paste0(root_url, endpoint)) |>
req_headers(Authorization = key) |>
req_perform()
print(resp_raw(test))
## HTTP/1.1 200 OK
## content-type: application/json
## X-Cloud-Trace-Context: 56fbd20bd5568e2528627f452c2d28bc;o=1
## Date: Wed, 31 Jan 2024 17:32:34 GMT
## Server: Google Frontend
## Content-Length: 35
## Via: 1.1 google
## Alt-Svc: h3=":443"; ma=2592000,h3-29=":443"; ma=2592000
##
## {"ClimateEngine":"Your key works!"}
## <httr2_response>
## GET https://api.climateengine.org//home/validate_key
## Status: 200 OK
## Content-Type: application/json
## Body: In memory (35 bytes)
# Clean up unneeded objects
rm(test, endpoint)
This section of the notebook runs using the Sandy Ecological Site Group, but you can run similar analysis for any shapefile. We will import the shapefile as an sf object and upload it to Google Earth Engine and then visualize the Sandy ESG on a leaflet map using the mapview package.
# Print the number of vertices in shapefile
# Because of the complexity of the shapefile, for this demo I uploaded it to Google Earth Engine
# For less complex shapefiles there are ways to pull the coordinates directly from the shapefile without Earth Engine
# Here is an example notebook: https://github.com/Google-Drought/SupportSiteTutorials/blob/a038ee1e69fff32008d8619c0acc6db082d5795d/Timeseries/OpenET_Example.Rmd
print(mapview::npts(sandy_esg))
## [1] 383672
# Visualize the allotments on a leaflet map using mapview package
mapview(sandy_esg)